In this section, we discuss in detail the steps required to analyze our empirical example and reproduce the results. You can run the code below to reproduce the results discussed in the empirical example section. The R Markdown file that contains all the code as well as the data set for the empirical example is available from a GitHub repository (https://github.com/quantPsych/mbco). When using the data set in this example, please cite the relevant study by MacKinnon, Valente, and Wurpts (2018).

Here are the packages required to reproduce the results in this document:

memory_df <- read.csv("memory_study.csv")
memory_df <-
memory_df %>% rename(repetition = R,
imagery = M,
recall = Y) %>%
mutate(x = factor(X, 0:1, labels = c("repetition", "imagery")))

We first gets a glimpse of the data set and it structure using the glimpse function from the dplyr package:

glimpse(memory_df)
## Observations: 369
## Variables: 7
## $ study      <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ X          <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
## $ repetition <int> 1, 1, 1, 3, 7, 2, 1, 7, 3, 6, 6, 9, 3, 1, 4, 4, 3, 7,…
## $ recall     <int> 9, 14, 17, 8, 10, 10, 8, 12, 12, 15, 7, 10, 8, 9, 15,…
## $ imagery    <int> 6, 8, 9, 7, 7, 6, 3, 6, 7, 7, 8, 8, 5, 5, 8, 9, 4, 7,…
## $ XM         <int> 6, 8, 9, 7, 7, 6, 3, 6, 7, 7, 8, 8, 5, 5, 8, 9, 4, 7,…
## $ x          <fct> imagery, imagery, imagery, imagery, imagery, imagery,…

Research Question 1: Does the Instruction to Create Mental Images of Words Increase Use of Mental Imagery that, then, Increases the Number of Words Recalled?

For this research question, we were interested in testing whether the indirect effect of Instruction of Recall through Imagery is different from zero. We used the single-mediator model (Figure 1) to answer this question. We followed the three steps outlined in the manuscript to implement the MBCO procedure. We also briefly introduced the key parts of OpenMx (Boker et al. 2011; Neale et al. 2016) code to implement the MBCO procedure in R.

Step 1

endVar <- c('imagery', 'recall') # name the endogenous variables
maniVar <-
  c('X', 'imagery', 'recall') # name the observed variables

single_med_full <- mxModel(
  "Single_med_full",
  type = "RAM",
  manifestVars = maniVar,
  #specify the manifest (observed) variables
  mxPath(
    from = "X",
    to = endVar,
    arrows = 1,
    free = TRUE,
    values = .2,
    labels = c("b1", "b3") # specify the path from X to the M (imagery) and Y (recall)
  ),
  mxPath(
    from = 'imagery',
    to = 'recall',
    arrows = 1,
    free = TRUE,
    values = .2,
    labels = "b2" # specify the path from M to Y
  ),
  mxPath(
    from = maniVar,
    arrows = 2,
    free = TRUE,
    values = .8,
    labels = c("s2x", "s2em", "s2ey") # specify (residual) variances for the observed variables
  ),
  mxPath(
    from = "one",
    to = endVar,
    arrows = 1,
    free = TRUE,
    values = .1, 
    labels = c("int1","int2") # specify the intercepts for the endogenous variables
  ),
  mxAlgebra(b1 * b2, name = "ind"), # define the indirect effect and call it ind
  mxData(observed = memory_df, type = "raw") # specify the data set for analysis
)
fit_single_med_full <- mxRun(single_med_full) # run the single mediator null model
## Warning: In model 'Single_med_full' Optimizer returned a non-zero status
## code 1. The final iterate satisfies the optimality conditions to the
## accuracy requested, but the sequence of iterates has not yet converged.
## Optimizer was terminated because no further improvement could be made in
## the merit function (Mx status GREEN).
stat_single_med_full <- summary(fit_single_med_full) # save the results
stat_single_med_full # print results of the full model
## Summary of Single_med_full 
##  
## The final iterate satisfies the optimality conditions to the accuracy requested, but the sequence of iterates has not yet converged. Optimizer was terminated because no further improvement could be made in the merit function (Mx status GREEN). 
##  
## free parameters:
##   name matrix row col   Estimate  Std.Error A
## 1   b1      A   2   1 3.64018505 0.24259998  
## 2   b3      A   3   1 0.04528978 0.38482959  
## 3   b2      A   3   2 0.58267974 0.06507646  
## 4  s2x      S   1   1 0.50948343 0.03750860  
## 5 s2em      S   2   2 5.42736402 0.39956681  
## 6 s2ey      S   3   3 8.48164638 0.62442551  
## 7 int1      M   1   2 3.80662759 0.17316342  
## 8 int2      M   1   3 8.74326503 0.32898072  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              8                   1099              4305.727
##    Saturated:              9                   1098                    NA
## Independence:              6                   1101                    NA
## Number of observations/statistics: 369/1107
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       2107.727               4321.727                 4322.127
## BIC:      -2190.239               4353.013                 4327.632
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2019-09-30 21:22:42 
## Wall clock time: 0.07884407 secs 
## optimizer:  NPSOL 
## OpenMx version number: 2.14.11.15 
## Need help?  See help(mxSummary)
plot(fit_single_med_full, means=FALSE, resid = "line", digits = 3)
## [1] "?plot.MxModel options: std, means, digits, strip_zero, file, splines=T/F, min=, max =, same = , fixed, resid= 'circle|line|none'"
DiagrammeR::grViz("Single_med_full.gv")

In the first two lines of the script above, before specifying the model parts with mxModel, we can simplify the process by grouping the names of variables as a vector that will be used in model specification. For example, we specified a vector of the names of observed variables and saved them as maniVar. Then, we specified a vector of the names of endogenous variables and saved them as endVar. Next, we specified the full single-mediator model.

The main command to specify an SEM is mxModel. The arguments provided to the mxModel function specified all the elements of the mediation model. For the single-mediator example, we specified the manifest (observed) and latent variables, the paths (regression coefficients), the indirect effect to be tested (or any function of model parameters), the constraint distinguishing the full and null models, and the variances and covariances among the variables.

The first argument to mxModel is single_med_full, which is a name we chose for the single-mediator model. The next argument, type=“RAM”, specifies that OpenMx uses the reticular action model (RAM; McArdle and McDonald (1984)), a symbolic algebraic notation to specify an SEM. In the argument manifestVars, we introduced the vector of the names of the observed (manifest) variables. The variable names must match the names in the data set memory_df, as shown earlier in the output from the glimpse function.

Next, we specify the paths between the variables using mxPath(). The function mxPath() corresponds to the graphical representation of paths in an SEM. For example, we use mxPath()to indicate a path (coefficient) corresponding to an arrow between the two variables specified in the arguments from (predictors) and to (response variables) in a single-mediator model. The argument arrows=1 indicates a unidirectional arrow that starts from the variable in the argument from and ends at the variable specified in the argument to; that is, a unidirectional arrow indicates a path coefficient between the two variables. The argument arrows=2 indicates a bidirectional arrow representing a covariance between the two variables. The argument free = TRUE indicates that the parameter is freely estimated; otherwise, free = FALSE indicates that the parameter is fixed at the values set by the argument values. If the parameter is freely estimated, the argument values would provide starting values. The argument labels provides labels for the coefficients. Because we specified more than one coefficient in the arguments toand from, we would provide more than one label corresponding to the stated order of the coefficients. For our example, b1 is the coefficient for \(X \rightarrow\) Imagery and b3 is the coefficient for \(X \rightarrow\) Recall.

We use mxAlgebra() to define the indirect effect or, in general, a function of model parameters. In general, a function may include mathematical operations +,-,*,/ (e.g., b1*b2/(b1*b2+b3)), exponential (e.g., exp()), and logarithms (e.g., log()). The first argument to mxAlgebra() is the product of two coefficients, b1 * b2, in which b1 and b2 had been defined in mxPtah(). The argument name = "ind" is used to name the indirect effect.

Next, we specified the data set for the model. The mxData identified the data set to be analyzed. The argument observed=memory_df specified the name of the data set in R. The second argument, type="raw", indicated that the data set was in the raw format, which meant that the data set included observations on the participants as opposed to being a summary statistic such as a covariance matrix.

Finally, we run the model using mxRun(), where the first argument is the name of the mxModel that is then saved as fit_single_med_full. After getting no warning about whether the model estimation and convergence criterion are satisfied, we use the function summary() to save or print the summary of the results. We save the summary of the results as stat_single_med_full and then print the summary. Below, we present the relevant part of the summary results.

Model Statistics:
             | Parameters | Degrees of Freedom | Fit (-2lnL units)
       Model:        8           1099             4305.727
   Saturated:        9           1098                NA
Independence:        6           1101                NA
Information Criteria:
     | df Penalty  | Parameters Penalty | Sample-Size Adjusted
AIC:     2107.727       4321.727             4322.127
BIC:    -2190.239       4353.013             4327.632

Below the title Model Statistics, the row that starts with Model, gives the pertinent information for the full model. The output shows that the full single-mediator model has eight free parameters, with \(df_{Full}=1099\) and \(\mathcal{D}_{Null}=\) 4305.726. Under Information Criteria and Parameters Penalty, we get the estimates of the information fit indices. For the full model, the information fit indices are AICFull= 4321.727 and BICFull= 4353.013. Alternatively, we can compute deviance, the AIC, and the BIC using the following functions:

-2*as.vector(logLik(fit_single_med_full))
## [1] 4305.727
AIC(fit_single_med_full)
## [1] 4321.727
BIC(fit_single_med_full)
## [1] 4353.013

Checking Model Assumptions

Before proceeding further, after fitting a mediation model, it is important to check the statistical assumptions about normality of the residuals and the presence of outliers (Cohen et al. 2003) . Because OpenMx does not offer checking residuals facilities, we fit the single-mediator model using two regression equations with the lm() function.

lm_image <- lm(imagery ~ x, data = memory_df)
lm_recall <- lm(recall ~ x + repetition, data = memory_df)
ggResidpanel::resid_panel(lm_image)

ggResidpanel::resid_panel(lm_recall)

outlierTest(lm_image)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##      rstudent unadjusted p-value Bonferroni p
## 285 -2.792633          0.0055027           NA
outlierTest(lm_recall)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##     rstudent unadjusted p-value Bonferroni p
## 248 2.848921          0.0046354           NA
influenceIndexPlot(lm_image)

influenceIndexPlot(lm_recall)

We checked for normality of the residuals using QQ plots which indicated that the normality assumption for the residuals was reasonable. We also checked for the outliers using the influence plots and \(t\)-tests (Fox 2016). The results showed there were no outliers.

Step 2

We ran the null (restricted) model, which is the single-mediator model in which the indirect effect of Instruction on Recall through Imagery is constrained to zero, \(\beta_1 \beta_2 =0\). The code for the null model shown below.

single_med_null <-
  mxModel(model = single_med_full,
          name = "single_med_null",
          mxConstraint(ind == 0, name = "b1b2_equals_0"))
fit_single_med_null <- mxRun(single_med_null)
stat_single_med_null <- summary(fit_single_med_null)
stat_single_med_null
## Summary of single_med_null 
##  
## free parameters:
##   name matrix row col     Estimate    Std.Error A
## 1   b1      A   2   1 9.522179e-23 1.038197e-23 !
## 2   b3      A   3   1 4.529096e-02 3.656067e-01  
## 3   b2      A   3   2 5.826781e-01 6.352900e-02  
## 4  s2x      S   1   1 5.094816e-01 3.749345e-02 !
## 5 s2em      S   2   2 8.738891e+00 6.423770e-01  
## 6 s2ey      S   3   3 8.481638e+00 6.234669e-01  
## 7 int1      M   1   2 5.661250e+00 1.538597e-01  
## 8 int2      M   1   3 8.743274e+00 3.282580e-01  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:              8                   1100              4481.493
##    Saturated:              9                   1099                    NA
## Independence:              6                   1102                    NA
## Number of observations/statistics: 369/1108
## 
## Constraint 'b1b2_equals_0' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       2281.493               4497.493                 4497.893
## BIC:      -2020.384               4528.779                 4503.398
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2019-09-30 21:22:45 
## Wall clock time: 0.04969096 secs 
## optimizer:  NPSOL 
## OpenMx version number: 2.14.11.15 
## Need help?  See help(mxSummary)

Note that instead of specifying all the parts of the null mediation model using MxModel in the above code, we modified the full model single_med_full, as specified in the argument model. Next, we specified the non-linear constraint for the indirect effect through mxConstraint. The first argument ind == 0 constrains the indirect effect defined in the mxAlgebra statement to zero. The argument name assigns a name to the constraint. Finally, we saved the null model to single_med_null. Next, we ran the null model and saved the results to fit_single_med_null. Relevant parts of the summary of the model results are shown below:

Model Statistics: 
               |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
       Model:              8                   1100              4481.493
   Saturated:              9                   1099                    NA
Independence:              6                   1102                    NA
Number of observations/statistics: 369/1108
Information Criteria: 
      |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       2281.493               4497.493                 4497.893
BIC:      -2020.384               4528.779                 4503.398

For the null model, \(df_\text{Null}= 1100\) and \(\mathcal{D}_\text{Null}=\) 4481.492, and the information fit indices were AICNull = 4497.493 and BICNull = 4528.779.

Step 3

We compared the full and null model both in terms of the LRTMBCO and the information fit indices to evaluate \(H_0: \beta_1 \beta_2 =0\). The LRTMBCO equals the difference between the deviance of the two models. We use the function code mxCompare to compute the LRTMBCO as follows:

lrt_rq1 <- mxCompare(fit_single_med_full, fit_single_med_null) 
lrt_rq1

The first row of the above output shows the results for the full model, which is the single-mediator model in Figure 1. The columns ep, minus2LL, df, and AIC show the number of estimated parameters, deviance, degrees of freedom, and AIC, respectively. The columns diffLL, diffdf, and p represent the difference in deviance (not the log-likelihoods), the difference in degrees of freedom, and p-value for two models being compared. The second row shows the results for the null model under the columns ep, minus2LL, df, and AIC. The results of comparing the null and full model are shown under the columns diffLL, diffdf, and p. The value for the test statistic LRTMBCO = 175.766 is located under the column diffLL. The degrees of freedom are dfLRT = 1 and p-value = 4.073738410^{-40}; these amounts are located under the columns diffdf and p, respectively.

Next, we used the following commands to compute the Monte Carlo CI for the indirect effect. First, we used the functions coef and vcov to extract the path coefficients and covariance matrix of the coefficients, respectively, from the full single-mediator model. Next, we used the ci function in the RMediation package (Tofighi and MacKinnon 2011, 2016) to compute the 95% Monte Carlo CI. The first argument mu to this function is a vector of the coefficient estimates, and the second argument Sigma is a covariance matrix of the coefficient estimates. The argument quant accepts a formula for the indirect effect that starts with the symbol “~”.

Mu <- coef(fit_single_med_full) # path coefficient estimates from the full single mediator model
Sigma <-
  vcov(fit_single_med_full) #covariance matrix of the parameter estimates from the full single mediator model

mc_ci1 <- RMediation::ci(mu = Mu,
                         Sigma = Sigma,
                         quant = ~ b1 * b2)
mc_ci <- t(round(unlist(mc_ci1[1:3]),3 ))
cat(" Monte Carlo CI for b1*b2:\n")
##  Monte Carlo CI for b1*b2:
knitr::kable(mc_ci)
2.5 % 97.5 % Estimate SE
1.6 2.682 2.121 0.276

We also recommend computing the difference in \(R^2\)s between the full and null model to examine change in the effect sizes that occurs as a result of the indirect effect through Imagery. As shown below, for Recall, \(R^2\) remained unchanged to four decimal places while for Imagery, \(\Delta R^2 = .3779\).

cat("R-Square for the full model\n")
## R-Square for the full model
rsq(model = fit_single_med_full, name = c("imagery","recall") )
##   imagery    recall 
## 0.3806254 0.2642803
cat("R-Square for the null model\n")
## R-Square for the null model
rsq(model = fit_single_med_null, name = c("imagery","recall") )
##     imagery      recall 
## 0.002711604 0.264281019
cat("Delta R-Square \n")
## Delta R-Square
round(rsq(model = fit_single_med_full, name = c("imagery","recall") )-rsq(model = fit_single_med_null, name = c("imagery","recall") ), 4)
## imagery  recall 
##  0.3779  0.0000

The first argument model to rsq specifies an OpenMx model and the second argument name specifies names of endogenous variables (i.e., variables in which a single headed arrow enters; namely those that are a function of another variable).

Research Question 2: Does the Instruction to Repeat Words Increase Use of Repetition to Memorize the Words that, in turn, Increases the Number of Words Recalled over and above Using Mental Imagery?

Step 1

We estimated the model with two parallel mediators in Figure 2, which is the full model for this research question. For this model, the two specific indirect effects associated with Imagery and Repetition were freely estimated. Below is the OpenMx code for the full parallel two mediator model:

endVar <- c('imagery', 'repetition', 'recall')
maniVar <- c('X', 'imagery', 'repetition', 'recall')
two_med_full <- mxModel(
  "two_med_full",
  type = "RAM",
  manifestVars = maniVar,
  mxPath(
    from = "X",
    to = endVar,
    arrows = 1,
    free = TRUE,
    values = .2,
    labels = c("b1", "b3", "b5")
  ),
  mxPath(
    from = 'repetition',
    to = 'recall',
    arrows = 1,
    free = TRUE,
    values = .2,
    labels = 'b4'
  ),
  mxPath(
    from = 'imagery',
    to = 'recall',
    arrows = 1,
    free = TRUE,
    values = .2,
    labels = "b2"
  ),
  mxPath(
    from = maniVar,
    arrows = 2,
    free = TRUE,
    values = .8,
    labels = c("s2x", "s2em1", "s2em2", "s2ey")
  ),
  mxPath(
    from = 'imagery',
    to = 'repetition',
    arrows = 2,
    free = TRUE,
    values = .2,
    labels = "cov_m1m2"
  ),
  mxPath(
    from = "one",
    to = endVar,
    arrows = 1,
    free = TRUE,
    values = .1,
    labels = c("int1", "int2", "int3")
  ),
  mxAlgebra(b1 * b2, name = "ind1"),
  mxAlgebra(b3 * b4, name = "ind2"),
  mxData(observed = memory_df, type = "raw")
)

fit_two_med_full <- mxTryHard(two_med_full) # run the full model
stat_two_med_full <- summary(fit_two_med_full) # save the results
stat_two_med_full # printing results
## Summary of two_med_full 
##  
## free parameters:
##        name matrix row col    Estimate  Std.Error A
## 1        b1      A   2   1  3.64017190 0.24259409  
## 2        b3      A   3   1 -3.80380675 0.21862053  
## 3        b5      A   4   1  0.10901983 0.44349780 !
## 4        b2      A   4   2  0.58769683 0.06734947  
## 5        b4      A   4   3  0.02155688 0.07473823  
## 6       s2x      S   1   1  0.50948566 0.03750880  
## 7     s2em1      S   2   2  5.42737309 0.39953786  
## 8  cov_m1m2      S   2   3 -1.26366712 0.26291464 !
## 9     s2em2      S   3   3  4.40768767 0.32447306  
## 10     s2ey      S   4   4  8.47970473 0.62426883 !
## 11     int1      M   1   2  3.80663024 0.17316070  
## 12     int2      M   1   3  8.01656645 0.15604890  
## 13     int3      M   1   4  8.55135498 0.74226839  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:             13                   1463              5874.685
##    Saturated:             14                   1462                    NA
## Independence:              8                   1468                    NA
## Number of observations/statistics: 369/1476
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       2948.685               5900.685                 5901.711
## BIC:      -2772.810               5951.526                 5910.281
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2019-09-30 21:22:47 
## Wall clock time: 0.02624488 secs 
## optimizer:  NPSOL 
## OpenMx version number: 2.14.11.15 
## Need help?  See help(mxSummary)
## Compute R2s for the endogenous variables
rsq(fit_two_med_full, c("imagery", "repetition", "recall"))
##    imagery repetition     recall 
##  0.3806244  0.4521506  0.2644487

Note that in above code, both endVar and maniVar now include the names of the mediators Imagery and Repetition. As a result, we specified two extra paths from Instruction to Repetition and from Repetition to Recall using mxPath. Because the two mediators can covary, we specified the covariance between the associated residuals and named it cov_m1m2. Finally, we specified two specific indirect effects through the two mediators using mxAlgebra and named them ind1 and ind2.

The above output shows that the full two-mediator model had 13 free parameters, with dfFull =1463 and \(\mathcal{D}\)Full = 5874.685. For the full model, the results were that AICFull = 5900.685 and BICFull = 5951.526; the effect sizes were R2Imagery =.38, R2Repetition =.45, and R2Recall =.26. We computed the specific indirect effects through Imagery and Repetition and the 95% CI for each specific indirect effect using the ci function in the RMediation package. The results are shown below:

Mu <- coef(fit_two_med_full) # path coefficient estimates
Sigma <-
  vcov(fit_two_med_full) #covariance matrix of the parameter estimates
## Monte Carlo CI for b1*b2:
mc_ci_b1b2 <- ci(mu = Mu,
             Sigma = Sigma,
             quant = ~ b1 * b2)
mc_ci_b1b2 <- t(round(unlist(mc_ci_b1b2)[1:3],3))
cat(" Monte Carlo CI for b1*b2:\n")
##  Monte Carlo CI for b1*b2:
knitr::kable(mc_ci_b1b2)
2.5 % 97.5 % Estimate
1.603 2.715 2.14 
## Monte Carlo CI for b4*b5:
mc_ci_b4b5 <- ci(mu = Mu,
             Sigma = Sigma,
             quant = ~ b3 * b4)
mc_ci_b4b5<- t(round(unlist(mc_ci_b4b5)[1:3],3))
cat(" Monte Carlo CI for b3*b4:\n")
##  Monte Carlo CI for b3*b4:
knitr::kable(mc_ci_b4b5)
2.5 % 97.5 % Estimate
-0.642 0.475 -0.082

Step 2

Next, we fit the null model that satisfies the null hypothesis: \(H_0: \beta_3 \, \beta_4 =0\). We fit the null model by fixing the mediation chain through Repetition to zero, \(\beta_3 \, \beta_4 =0\). In OpenMx, we specified the null model by adding the constraint mxConstraint(ind_repetition == 0, name = "b4b5_equals_0") to the full model in Step 1, ran it, and then saved the results to fit_two_med_null. The summary results as well as the effect sizes for the endogenous variables for the null two-mediator model are shown below:

two_med_null <- mxModel(model = two_med_full,
                        name = "two_med_null",
                        mxConstraint(ind2 == 0, name = "b3b4_equals_0"))
fit_two_med_null <- mxRun(two_med_null)
stat_two_med_null <- summary(fit_two_med_null) # saving the results
stat_two_med_null # printing results 
## Summary of two_med_null 
##  
## free parameters:
##        name matrix row col      Estimate    Std.Error A
## 1        b1      A   2   1  3.640174e+00 2.416781e-01  
## 2        b3      A   3   1 -3.803803e+00 2.181357e-01  
## 3        b5      A   4   1  4.528902e-02 3.630610e-01  
## 4        b2      A   4   2  5.826797e-01 6.331447e-02  
## 5        b4      A   4   3  5.129810e-23 2.941779e-24 !
## 6       s2x      S   1   1  5.094832e-01 3.749198e-02  
## 7     s2em1      S   2   2  5.427373e+00 3.966568e-01  
## 8  cov_m1m2      S   2   3 -1.263665e+00 2.574043e-01  
## 9     s2em2      S   3   3  4.407693e+00 3.221346e-01  
## 10     s2ey      S   4   4  8.481652e+00 6.233655e-01  
## 11     int1      M   1   2  3.806630e+00 1.726238e-01  
## 12     int2      M   1   3  8.016564e+00 1.558209e-01  
## 13     int3      M   1   4  8.743264e+00 3.281558e-01  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:             13                   1464              5874.768
##    Saturated:             14                   1463                    NA
## Independence:              8                   1469                    NA
## Number of observations/statistics: 369/1477
## 
## Constraint 'b3b4_equals_0' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       2946.768               5900.768                 5901.794
## BIC:      -2778.638               5951.609                 5910.364
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2019-09-30 21:22:51 
## Wall clock time: 0.03534794 secs 
## optimizer:  NPSOL 
## OpenMx version number: 2.14.11.15 
## Need help?  See help(mxSummary)
## Compute R2s for the endogenous variables
rsq(fit_two_med_null, c("imagery", "repetition", "recall"))
##    imagery repetition     recall 
##  0.3806244  0.4521500  0.2642798

Step 3

Finally, we compared the two models in terms of the LRTMBCO:

mxCompare(fit_two_med_full, fit_two_med_null)
# Compute differences in R2 for the endogenous variables
round(rsq(fit_two_med_full, c("imagery", "repetition", "recall"))- rsq(fit_two_med_null, c("imagery", "repetition", "recall")), 4)
##    imagery repetition     recall 
##      0e+00      0e+00      2e-04

The results showed that LRTMBCO = 0.083, dfLRT =1, and p = .773. The specific indirect effect through Repetition was, therefore, not different from zero, -0.08 (SE = 0.29), 95% Monte Carlo CI = [-0.64, 0.48]. Further, the R2 for Imagery, Repetition, and Recall remained unchanged to three decimal places, and the information fit indices between the two models were roughly the same. These results indicate that the specific indirect effect through Repetition above and beyond the specific indirect effect through Imagery does not appear to be different from zero.

Research Question 3: Is the Indirect Effect of Instruction to Create Mental Images on the Number of Words Recalled through the Use of Mental Imagery Greater than the Indirect Effect of Instruction to Repeat Words on the Number of Words Recalled Through the Use of Repetition?

For the third research question, we were interested in comparing the sizes of the two specific indirect effects: the indirect effect of Instruction on Recall through Imagery (i.e., \(β_1 β_2\)) and the indirect effect of Instruction on Recall through Repetition (i.e., \(β_4 β_5\)).

Step 1

The full model for this research question was the same as the full model in Research Question 2. Thus, we used the results (i.e., the vector of the coefficient estimates, the covariance matrix of the coefficient estimate, R2, indirect effects estimates, AIC, and BIC) of the full parallel two-mediator model from Research Question 2.

Step 2

The null hypothesis for this research question is \(H_0: \beta_1 \beta_2 = \beta_3 \beta_4\). To test this hypothesis, the null model is a two-mediator model in which we constrained the two specific indirect effects to be equal. To specify the null model for this research question in OpenMx, we added the following argument to the mxModel function for full model in Research Question 2: mxConstraint(ind_imagery==ind_repetition,name = "ind1_eq_ind2"). We then ran the model and saved the results as fit_two_med_contrast.

## two-mediator model with the contrast fixed at zero
two_med_contrast <- mxModel(model = two_med_full,
                            name = "two_med_contrast",
                            mxConstraint(ind1 == ind2, name = "ind1_eq_ind2"))
fit_two_med_contrast <- mxRun(two_med_contrast) # fitting the model
stat_two_med_contrast <-
  summary(fit_two_med_contrast) #save the results
stat_two_med_contrast ## print the summary results
## Summary of two_med_contrast 
##  
## free parameters:
##        name matrix row col   Estimate  Std.Error A
## 1        b1      A   2   1  3.4619985 0.23928899 !
## 2        b3      A   3   1 -3.9128023 0.21720490 !
## 3        b5      A   4   1 -0.1882543 0.40821963  
## 4        b2      A   4   2  0.3360970 0.04573403 !
## 5        b4      A   4   3 -0.2973744 0.04132634 !
## 6       s2x      S   1   1  0.5094833 0.03748574  
## 7     s2em1      S   2   2  5.4353060 0.39613042  
## 8  cov_m1m2      S   2   3 -1.2588089 0.25480530  
## 9     s2em2      S   3   3  4.4106633 0.32120360  
## 10     s2ey      S   4   4  9.0688068 0.66792636  
## 11     int1      M   1   2  3.8974080 0.17176855  
## 12     int2      M   1   3  8.0720977 0.15548889  
## 13     int3      M   1   4 12.0658321 0.29118476  
## 
## Model Statistics: 
##                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
##        Model:             13                   1464              5900.514
##    Saturated:             14                   1463                    NA
## Independence:              8                   1469                    NA
## Number of observations/statistics: 369/1477
## 
## Constraint 'ind1_eq_ind2' contributes 1 observed statistic. 
## 
## Information Criteria: 
##       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
## AIC:       2972.514               5926.514                 5927.539
## BIC:      -2752.893               5977.354                 5936.109
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2019-09-30 21:22:51 
## Wall clock time: 0.03577805 secs 
## optimizer:  NPSOL 
## OpenMx version number: 2.14.11.15 
## Need help?  See help(mxSummary)
# Compute R2s for the null model
rsq(fit_two_med_contrast, c("imagery", "repetition", "recall"))
##    imagery repetition     recall 
##  0.3797190  0.4517808  0.2133485

The above output shows that the null model had dfNull =1464 and DNull =5900.514. The information fit indices for the null model were AICNull = 5926.514 and BICNull = 5977.354. The effect sizes were R2Imagery =.38, R2Repetition=.45, and R2Recall =.21. # Step 3

We compared the full parallel two-mediator model in Step 1 and the null parallel two-mediator model in Step 2. We computed the LRTMBCO using the function mxCompare(fit_two_med_full,fit_two_med_contrast):

mxCompare(fit_two_med_full,fit_two_med_contrast) 
# Compute differences in R2s
round(rsq(fit_two_med_full, c("imagery", "repetition", "recall"))-rsq(fit_two_med_contrast, c("imagery", "repetition", "recall")), digits = 4)
##    imagery repetition     recall 
##     0.0009     0.0004     0.0511

We also computed the 95% Monte Carlo CI for the contrast of the two indirect effects using the ci function in the RMediation package, where the arguments mu and Sigma were the vector of the coefficient estimates and the covariance matrix of the coefficient estimates, respectively, obtained from the full parallel two-mediator model in Step 1:

Mu <- coef(fit_two_med_full) # path coefficient estimates
Sigma <-
  vcov(fit_two_med_full) #covariance matrix of the parameter estimates
mc_ci_con <- ci(mu = Mu,
                Sigma = Sigma,
                quant = ~ b1 * b2 - b3 * b4)
mc_ci_con <- t(round(unlist(mc_ci_con)[1:3],3))

cat(" Monte Carlo CI for Contrast:\n")
##  Monte Carlo CI for Contrast:
knitr::kable(mc_ci_con)
2.5 % 97.5 % Estimate
1.364 3.11 2.222

The results of the MBCO procedure showed that LRTMBCO = 25.828, dfLRT =1, and p=3.731857E-07. These outcomes indicate that the indirect effect through Imagery appeared to be larger than the indirect effect through Repetition by 2.222 (SE= 0.445) words, 95% Monte Carlo CI = [1.364, 3.11]. Comparing the R2s for the endogenous variables obtained from Step 1 and 2 for Imagery and Repetition, R2s remained unchanged to three decimal places while for Recall \(\Delta R^2\)=.05. Comparing the information fit indices of the null model to the full model also supports the conclusion that the full model fit the data better than did the null model. We next describe simulation studies that compare the statistical properties of LRTMBCO with the commonly used methods of testing indirect effects.

References

Boker, Steven, Michael C. Neale, Hermine Maes, Michael Wilde, Michael Spiegel, Timothy Brick, Jeffrey Spies, et al. 2011. “OpenMx: An Open Source Extended Structural Equation Modeling Framework.” Psychometrika 76 (2): 306–17. https://doi.org/10.1007/s11336-010-9200-6.

Cohen, Jacob, Patricia Cohen, Stephen G. West, and Leona S. Aiken. 2003. Applied Multiple Regression/Correlation Analysis for the Behavioral Sciences. 3rd ed. Mahwah, NJ: Erlbaum.

Fox, John. 2016. Applied Regression Analysis and Generalized Linear Models. 3rd ed. Los Angeles, CA: SAGE.

MacKinnon, David P., Matthew J. Valente, and Ingrid C. Wurpts. 2018. “Benchmark Validation of Statistical Models: Application to Mediation Analysis of Imagery and Memory.” Psychological Methods 23: 654–71. https://doi.org/10.1037/met0000174.

McArdle, J. Jack, and Roderick P. McDonald. 1984. “Some Algebraic Properties of the Reticular Action Model for Moment Structures.” British Journal of Mathematical and Statistical Psychology 37 (2): 234–51. https://doi.org/10.1111/j.2044-8317.1984.tb00802.x.

Neale, Michael C., Michael D. Hunter, Joshua N. Pritikin, Mahsa Zahery, Timothy R. Brick, Robert M. Kirkpatrick, Ryne Estabrook, Timothy C. Bates, Hermine H. Maes, and Steven M. Boker. 2016. “OpenMx 2.0: Extended Structural Equation and Statistical Modeling.” Psychometrika 81: 535–49. https://doi.org/10.1007/s11336-014-9435-8.

Tofighi, Davood, and David P. MacKinnon. 2011. “RMediation: An R Package for Mediation Analysis Confidence Intervals.” Behavior Research Methods 43: 692–700. https://doi.org/10.3758/s13428-011-0076-x.

———. 2016. “Monte Carlo Confidence Intervals for Complex Functions of Indirect Effects.” Structural Equation Modeling: A Multidisciplinary Journal 23: 194–205. https://doi.org/10.1080/10705511.2015.1057284.